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Abstract 

An algorithm for simulation of quantum many-body dynamics having su(2) spectrum-generating 
algebra is developed. The algorithm is based on the idea of dynamical coarse-graining. The 
original unitary dynamics of the target observables — the elements of the spectrum-generating 
algebra — is simulated by a surrogate open-system dynamics, which can be interpreted as weak 
measurement of the target observables, performed on the evolving system. The open-system state 
can be represented by a mixture of pure states, localized in the phase-space. The localization 
reduces the scaling of the computational resources with the Hilbert space dimension n by factor 

3 /2 

compared to conventional sparse- matrix methods. The guidelines for the choice of parameters 
for the simulation are presented and the scaling of the computational resources with the Hilbert 
space dimension of the system is estimated. The algorithm is applied to the simulation of the 
dynamics of systems of 2 • 10 4 and 2 • 10 6 cold atoms in the double-well trap, described by the 
two-sites Bose-Hubbard model. 

PACS numbers: 03.67.Mn,03.67.-a, 03.65.Ud, 03.65 Yz 
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I. INTRODUCTION 



Simulation of quantum many-body dynamics is a challenging task. The computational 
resources of solving the Schrodinger equation increase with the size of the system's Hilbert 

n 

space n, which is very large in a typical many-body Computational methods based 
on diagonalization scale as the the cube of the Hilbert space dimension. This led to the 
development of propagation techniques which scale as 0(nhxn) ( semi-linear ly) with the 
Hilbert space dimension for an elementary propagation time-step. The number of time- 
steps required for the simulation on a fixed time interval scales with the energy uncertainty 
of the evolving state, which adds generically at least another factor n to the computational 
cost. As a result, the overall scaling of the computational resources for a fixed time interval 
becomes O (n 2 In n) . To achieve such a scaling different sparse-matrix techniques are required 
such as FFT |3| or linear scaling approaches in solid state |4|. Currently such methods are 
able to address quantum dynamics in Hilbert spaces up to dimension of ~ 10 8 . Eventually 
due to the unfavorable scaling the semi-linear techniques become unpractical. 

As a consequence, other strategies have been developed to meet the challenge of the 
many-body quantum simulations. A popular approach is based on extensions of a mean-field 
approximation and is restricted to systems with a limited scope of quantum correlations. An 
example is the powerful multi configuration time dependent Hartree (MCDTH) algorithm, 
where part of the correlations are deliberately added back [5|. As the quantum correlations 
develop in the course of the evolution the scaling of the method becomes prohibitive. 

The present paper describes the first implementation of an alternative approach to many- 
body quantum simulations [2], which focuses on the dynamics of a small subset of observ- 
ables. The motivation for the approach is to find a method of quantum-dynamical simu- 
lations, scaling with the size of the subset of target observables but not with the Hilbert 
space dimension of the system. Simulations with such scaling properties are defined in the 
present work as efficient. The basic idea of the method proposed in Ref. is to simulate 
the unitary dynamics of a subset of observables by a surrogate open-system dynamics. The 
target observables in the proposed scheme are elements of the spectrum-generating algebra 
(SGA) of the quantum system j(| and the surrogate open-system dynamics corresponds to 
weak measurement [7| of the observables performed on the evolving system. 

The present study explores basic steps required to generate an operational algorithm and 
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the scaling of the computational resources with the Hilbert space dimension for simulation 
of quantum many-body dynamics having su(2) SGA. The expositions is carried out through 



the study o 



we. 



1 trap 



a specific example, the simulation of the dynamics of N cold atoms in a double- 



id, 
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131 ] . The system is modeled by the two-site Bose-Hubbard model 



14] with the Hamiltonian: 



H 



U 



~ A Yl $+1^ + + 2iV S 



(1) 

i=l,2 " ' i=l,2 

where Obi ) is the creation (annihilation) operator of a particle in the i-th well. A describes 
the hopping rate and U scales the two-body interaction. The transformation to the operators 



J, 



— (a{&2 + ai,ai y 



Jz = ^(^l^ 1 ~ ^2^2) 

leads to the following Lie-algebraic form of the Hamiltonian 

U 



(2) 



H 



(3) 



where u> = 2 A. The set of single-particle operators {J x , J y , J z } is closed under the commu- 
tation relation and spans the su(2) SGA of the system. Examples of many-body quantum 
systems with su(2) SGA include other systems such as the Lipkin-Meshkov-Glick model of 
inte^ingf— Q. The Hilbert s pac e of the system of iV atoms carries the n = N + 1 
dimensional irreducible representation [16| of the the algebra, corresponding to the pseudo- 
spin quantum number j = N/2. Observables J x , J y , J z are the target of the dynamical 
simulation, representing coherence, current and population difference between the wells. 

The corresponding surrogate open-system dynamics is driven by the Liouville-Von Neu- 
mann equation of motion: 



P 



H,p 



7 



(4) 



The simulation of the open-system evolution (J3J) is performed by averaging over solutions 
of the stochastic nonlinear Schrddinger equation (sNLSE) 17|, [l8|, ll9| . The sNLSE governs 



a single realization of the process of weak measurement of the target observables ([2]) per- 
formed on the system with the Hamiltonian ([3]). Parameter 7 in Eq.(j4]), corresponding to 



the strength of measurement, is decreased in the course of the simulation until the open 
dynamics (j3J) of the target observables converges to the original unitary dynamics, driven 
by the Hamiltonian Q. It is shown that the convergence is achieved at the magnitude of 7 
which is still sufficiently large to drastically reduce of the computational complexity of the 
sNLSE, compared to the original Schrodinger equation. The reduction of the computational 
complexity is achieved by virtue of the localization of solutions of the sNLSE in the associ- 
ated phase space, induced by the measurement 0, Q] . On the level of the density matrix 
the localization leads to a coarse-graining of the phase-space representation of the state, 
erasing the fine structure but leaving unaffected the dynamics of the target observables. 

Localization of the sNLSE solution leads to compression of the size of the time-dependent 
computational basis of the simulation to M < n = N+l states. In the convergent simulation 
on the time-interval of order Inn, M < Inn. As a consequence, the memory resources and 
the cost of an elementary time-step of the computation scale as In 2 n. The number of 
time-steps in the simulation on the time-interval of order Inn is 0(y/nlnn) for a localized 
state evolution. Averaging over realizations of the weak measurement adds a factor of 
~ 10 3 — 10 4 . This factor is independent on n. Therefore, theoretically, the overall scaling of 
the computational resources is 0(^Jn\v? n). Simulations using conventional sparse-matrix 
techniques scale as 0(n 2 In 2 n). Therefore, even though the algorithm is not efficient in the 
sense of our definition, the scaling of the computational resources is better by the factor of 
n 3 ' 2 (lnn) _1 than in conventional sparse-matrix approaches. 

An additional important feature of the algorithm is that simulations of the sNLSE for 
individual measurement realizations are independent and can be performed in parallel. This 
parallelism is uncommon in conventional algorithms for quantum dynamical simulations. 

Section Hit A) summarizes the main ideas of Refs. 2], [20] and motivates the choice of the 
computational basis for simulation of the target observables. Section [TT^B) develops a scheme 
for updating the computational basis using a global unitary transformation generated by 
the SGA . 

Section II II I focuses on the application of the algorithm to dynamics of a large number 
N of cold atoms in a double-well trap. The choice of parameters of the simulation and the 
scaling of the computational resources with the Hilbert space dimension are analyzed and 
numerical results of the simulation for N = 2 ■ 10 4 and for N = 2 ■ 10 4 atoms are presented. 
Section IIVI summarizes the results. 
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II. THE ALGORITHM 



A. General structure 

The algorithm proposed in Ref. 2J, is based on the following sequence of steps: 

1. Construction of the time-dependent computational basis of generalized coherent states 
(GCS), associated with the SGA of the system. 

2. The choice of elements of SGA of the system as the target observables for the simula- 
tion. 

3. Numerical solution of the sNLSE, corresponding to a single realization of weak mea- 
surement of the target observables, performed on the evolving system. 

4. Averaging the expectation values of the target observables over a large number of 
realizations. 

5. Convergence of the open-system dynamics of the target observables to their origi- 
nal unitary evolution, achieved by increasing the size of the computational basis and 
decreasing the strength of the measurement. 

1. The choice of the time- dependent basis 

As a guideline for the development of the algorithm the following definition of efficient 
simulation is adopted: 

Definition: The simulation of dynamics of a Lie-algebra of observables of the system is 
efficient if the necessary memory and the CPU resources for the computation do not depend 
on the Hilbert space (irreducible) representation of the algebra. 

In short, the simulation is called efficient if it is based on group-theoretical calculations. 
An example of efficient simulation is solution of the Heisenberg equations of motion for a 
Lie-algebra g of observables, when the Hamiltonian belongs to the algebra. In this case 
the simulation is reduced to numerical solution of dim {g} linear equations of motion for the 



expectation values of the elements o 



g. Another example is simulation of quantum dynamics 



in the mean-field approximation 21[ , when the computational complexity is independent on 



the Hilbert space dimension asymptotically. 



The starting point for the algorithm is the idea that reduction of the computational 
complexity can be achieved by using time-dependent basis of a small number M of states. 
Let ipiit) be the time-dependent computational basis. The state vector of the system can 
be represented in the time-dependent basis as 

M M 

!*(*)> = E C *W W*)> = X>(*)fy(*) M°)> , (5) 



i=l i=\ 



where M < n and Uj(i) is a non unique time-dependent unitary transformation of the 
reference state ipi(0). The simulation involves calculation of matrix elements of the Hamil- 
tonian in the time-dependent basis. An efficient computation of the matrix elements of the 
Hamiltonian \ip(t)i H ip(t)j) is possible only if both the transformation 



h -> ti^mjjit). (6) 

and the computation of matrix elements 

Vi(0) H ^(0)\ (7) 



can be performed efficiently. 

Let operators j span the SGA of the system. Then the Hamiltonian can be repre- 
sented as a polynomial in Xj, as in Eq.([3]). If the unitary transformation Uj(t) is generated 
by an element of the SGA, the transformation (jSJ) can be performed group-theoretically 16], 
i.e., efficiently. It follows that if each ipi{t) in Eq.Q can be represented as an orbit of the 
reference state ^(0) under the action of the SGA: 

A(t) = e«Si«iW^^.( ), (8) 
X, G SGA, (9) 

the transformation ([6]) can be performed efficiently. Eq. ([8|) defines a generalized coherent 



16 



22 1 . The choice of 



state (GCS), associated with the SGA and the reference state ^(0) 
the time-dependent computational basis of GCS, associated with the SGA of the system, 
guarantees that the transformation ([6]) can be performed efficiently. 

Complexity of computing matrix elements (JTj) generally depends on the Hilbert space 
representation of the SGA. But if the reference state ^(0) is a maximal (minimal) weight 
state of the representation, the computation can be performed group-theoretically. These 
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considerations suggest the choice of the basis of the GCS, associated with the SGA of the 
system and a maximal (minimal) weight states of the representation, as the time-dependent 
computational basis of the computation. 



2. The choice of the target observables for the simulation 

The choice of computational basis of GCS, associated with the SGA of the system and 
a maximal (minimal) weight state of the representation, allows computing matrix elements 
ip(t)jj of an observable X efficiently if X is a polynomial in Xj. Such observables 
can be used as target observables for the efficient simulation. Elements Xj of the SGA are 
the simplest choice of the target observables. 



3. Solving the stochastic Nonlinear Schrddinger Equation 

The size M of the computational basis in an efficient simulation must be independent 
of the Hilbert space dimension. Generically a unitarily evolving system will evolve into a 
superposition of M ~ n GCS. Therefore, generically, a unitary dynamics cannot be simu- 
lated efficiently. The open-system dynamics, corresponding to the weak measurement of the 
SGA operators of the evolving system, is simulated by averaging over stochastic pure-state 

ach realization is governed by the stochastic Nonlinear 



realizations of the measurement. 
Schrddinger equation (sNLSE)^. 



191 ] . Solution of the sNLSE, can be expressed in the 



form Eq.Q with M <C n, provided the measurement is sufficiently strong. 



4- Averaging 

Expectation values of the SGA observables are calculated in each stochastic realization 
of the sNLSE. Averaging over the realizations recovers the open-system dynamics of the 
observables. The number of realizations needed to compute the averages to a given rela- 
tive error is independent of the Hilbert space dimension. Therefore, the averaging can be 
performed efficiently. 



7 



5. Convergence 



Numerical solution of the sNLSE converges as the size M of the computational basis 
is increased. The open-system dynamics of the SGA observables converges to the original 
unitary dynamics as the strength 7 of the measurement is decreased. The smaller 7 the larger 
M is needed for the convergent simulation. If a certain generic condition on the SGA and its 
Hilbert space representation is satisfied |2|] , the open system follows the evolution on widely 
separated time-scales. Due to this time-scales separation the open-system dynamics of the 
SGA observables converges to the unitary evolution at the value of 7 which is sufficiently 
large to ensure that M <C n in the expansion (J5J) of the convergent pure-state solution of 
the corresponding sNLSE. 

B. A linear implementation 

The algorithm outlined in the previous subsection can be implemented in a variety of 
ways. This is in part due to the non-uniqueness of the transformations Uj(t) of the com- 
putational basis elements in Eq.(jSJ). In the present implementation we adopt the choice 
Ui(t) = Uj(t) = U(t), termed the linear implementation. In the linear implementation the 
pure-state solution of the sNLSE is represented by 



where ipi(0) is a fixed set of the GCS, associated with the SGA and a maximal (minimal) 
weight state of the Hilbert space representation of the algebra. 

The advantage of this choice is that updating of the computational basis is performed 
by a linear (unitary) transformation, leading to high numerical stability. The disadvantage 
is a limited flexibility of the computational basis evolution, which may result in a higher 
computational cost. 

An infinitesimal time step of the simulation can be represented by a superposition of the 
nonlinear transformation Ti of the state represented in the fixed instantaneous basis, i.e., 
the transformation of the vector Cj in Eq.f fTOl . and a subsequent unitary transformation T2 
of the basis itself. The transformation Ti corresponds to the evolution of the state, driven 
by the sNLSE. The transformation T 2 corresponds to the updating of the computational 



M 




(10) 



i=i 



s 



basis and ensures that the evolving state stays in the Hilbert subspace spanned by the 
instantaneous basis. Connection of the infinitesimal unitary transformation T2 to U(t) in 
Eq. lflUI) is established by the following relation: 

T 2 = l-\]{t)'dt. (11) 

The transformations T\ and T 2 can be given a simple geometric interpretation in the 
special case of the su(2) SGA spanned by the operators J x , J y and J z , satisfying the com- 
mutation relations 

Jj; Jj = i^ijk^k- (12) 

The GCS of the su(2), termed spin-coherent states, can be represented by points in the 



phase space, associated with the algebra [2J, |25|. The phase-space is a two-dimensional 



sphere and the GCS can be parametrized by spherical coordinates. The 2j + 1 dimensional 
Hilbert space, carrying an irreducible representation of the algebra, is spanned by a linearly 
independent basis of GCS, represented by 2j + 1 points on the sphere. A localized state can 
be expanded in a small fraction of the total basis, represented by points in the interior of a 
relatively small area of the phase-space. 

A single realization of the weak measurement of operators J x , J y and J z of a quantum 
system, evolving under the action of a Hamiltonian H is described by the following sNLSE: 

d = -z'H dt 

i=l ^ i=l ^ 

where (^ij are the instantaneous expectation values of the single-particle operators and 
the Wiener fluctuation terms d& satisfy 

= 0, [d&dtj] = 2-ySijdt, (14) 

where [ ] denotes statistical averaging. Parameter 7 is termed the strength of measurement 
in what follows. 

Hamiltonian H nonlinear in J x , J y and J z generates derealization of an initially localized 
state in the phase-space. The non-unitary stochastic contribution of the evolution due to the 
measurement, represented by the last two terms in Eq. ffl3|) . leads to localization of the state 



[20( | . A single realization of the open system evolution can be visualized in the phase-space 
ClS db stochastic trajectory of finite width. 

The infinitesimal nonlinear transformation Ti 

% = i + l-mdt (j^) dt + 2 ( 5i - (^} J ' ( 15 ) 

includes two parts: the original Hamiltonian part and the non-unitary part, which contains 
a stochastic and a deterministic element. Both elements of the non- unitary part depend on 
the instantaneous state of the system. 

The transformation T\ generates a drift of the state in the phase-space. The computa- 
tional basis must be updated to compensate for the drift. The GCS basis is updated at each 
time step by performing the unitary transformation T 2 , Eq. (II II) . Practically, a different 
but physically equivalent transformation has been found to be more convenient. Instead 
of updating the basis to follow the evolving state we update the state to match the fixed 
basis. Therefore, the inverse unitary transformation T 2 is performed on the state itself and 
the transformation T 2 is applied to the Hamiltonian and other observables. The unitary 
transformation T 2 is generated by the su(2) elements J x , J y and J z and rotates the state to 
the position localized in the fixed basis of the GCS. 

For example, let the initial state be a spin-coherent state, represented by the south pole of 
the phase space. The expectation values of the operators J x and J z vanish in this state. An 
infinitesimal transformation Tj takes the state into a new position having finite expectation 
values of J x and J z . The purpose of the unitary transformation T], is to rotate the state back 
to the state where the expectation values of J x and J z vanish. This rotation is performed 
by the 577(2) group transformation and is equivalent to rotating the density operator in the 
Hilbert-Schmidt space to the state where the projection of the density operator on the SGA 
is spanned by the J z operator. This procedure is analogous to the Schmidt-decomposition 



23]. 



of the state in the generalized entanglement setting 

The back-rotation of the operators by the transformation T 2 can be performed analyti- 
cally, knowing how the su(2) algebra elements transform under the 577(2) group itself. The 
superposition of the two transformations T 2 o Ti of the state results in a new state, localized 
about the south pole. 
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III. APPLICATION TO THE TWO-SITE BOSE-HUBBARD MODEL 



A. Introduction 

The gas of cold bosonic atoms in a double-well trap is described by the two-site Bose- 
Hubbard model, Eq.flT]). The Hamiltonian of the system contains two terms: the one-body 
term responsible for hopping of the atoms from one well to another and the two-body or the 
interaction term responsible for the on-site attraction or repulsion of the atoms. Here we 
consider repulsive interaction, measured by the coupling strength 211/ N, U > 0, where N is 
the total number of particles in the trap. The hopping rate is determined by the coefficient 
A of the one-body part. 

The SGA of the system is the su(2) algebra of the single-particle operators (J2J). The 
system of N particles correspond to the N + 1-dimensional irreducible representation of the 
algebra with the principle pseudo-spin quantum number j = N/2. Condensed state of the 



system is a spin coherent state 12 



13|. 



The character of dynamics depends on a single parameter U/u, where u = 2 A. At 
N ^> 1 the weak interaction regime U/u < 1, corresponds to dynamics preserving coherence 
of the atomic condensate for times of order \fN [26(. High coherence corresponds to a 
weak derealization of an initial spin-coherent state. The corresponding dynamics is very 
accurately described by the mean-field solution of the Schrodinger equation [Cf. Fig. QQ, 
assuming a spin-coherent form of the evolving state 211 ] . 

When U/u > 1 the phase-space picture of the mean-field dynamics changes: an unstable 
fixed point appears and the associated separatrix [3] . The mean-field solution breaks down 
in the vicinity of the separatrix. As a consequence, an initial spin-coherent state prepared in 
the vicinity of the separatrix is expected to evolve into a delocalized state. The corresponding 
dynamics is characterized by the exponential loss of coherence {27] as measured by the purity 
of the single-particle density operator. The derealization corresponds to depletion of the 
condensate. This parametric regime poses the real challenge for the simulation and is chosen 
to test the algorithm. 



The loss of coherence can be measured by the generalized purity 



231 ] of the state with 



respect to the su(2) algebra of the single-particle observables J x , J y and J 2 . The generalized 
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FIG. 1: (Color online) Mean- field (dashed lines) vs. exact (solid lines) solutions of the Schrodinger 
equation, corresponding to the Hamiltonian ([3]). Expectation values of the target observables 
= jf Jx (blue), J' y = jf Jy (red) and J' 2 = -^J z (green), Eq.(|2]), are plotted as functions of time 
for N = 128, N = 256 and N = 512 atoms. The generalized purity (black), Eq. (|16p . is plotted in 
N = 512 case. The system is prepared in a spin-coherent state corresponding to the condensate 
placed in the left well of the trap. The value of the on-site interaction is chosen U/u = 1. Time 
is measured in units of uj~ 1 . As the number of particles grows the exact solution approaches the 
mean-field solution. 



purity of the state ip is defined by 



P«MM = j!j5 ^ W- ( 16 ) 

i=x,y,z 

The generalized purity equals unity if ip is a spin-coherent state and is less than unity 
otherwise. Since the mean-field state is spin-coherent by construction, the mean-field ap- 
proximation breaks down if the generalized purity decreases. 

Figs. [T] and [2] show dynamics of the system of iV = 128, iV = 256 and iV = 512 atoms 
corresponding to U/u — 1 and U/u = 2 respectively, obtained by solving the Schrodinger 
equation using the Chebyshev method 28[. For initial state of the system we choose a 
condensate of atoms prepared in the left well. It is a spin-coherent eigenstate of the operator 
J z . The expectation value of J 2 measures the difference of the right- and left- well populations 



12 



time in units of co 



-1 



FIG. 2: (Color online) The same as in Fig(TJ but for U/u = 2. The mean-field approximation 
breaks down. 

of the atoms. The mean- field solution works well for the weak interaction regime U/u = 1, 
where the generalized purity of the evolving remains close to unity. In the case of U/u — 2 
the same initial state follows dynamics in the vicinity of the separatrix in the mean-field 
picture As a consequence, the generalized purity decreases on the time-scale u; _1 miV 
following the depletion of the condensate and the mean-field solution fails to reproduce 
correctly the character of the dynamics. 

B. Application to a system of large number of atoms 

1. Equations of motion 

The numerical solution of the sNLSE ({TBI) proceeds as a sequence of infinitesimal trans- 
formations: 





X ~+ T^XT; 



2- 



(17) 
(18) 
(19) 
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where X stands for all relevant operators, i.e., the Hamiltonian and the target observables 
P). The state ip in Eqs.fllZH and (ffBD is S iven b Y 

M 



Mt)) = E c ^) w°)> 



(20) 



i=l 



where ^(0) is the /ixerf computational basis of spin-coherent states. From Eqs. ffTTD . (fT5j) 
and (ED we obtain 



At 



Cjfc ^ 

fc=i 

k=l I 1=1 ^ ' 1=1 ^ 



(21) 



where (X) ifc = (^j(O) X ^ fc (0)> and 



C iJfc = <Vi(0)|^(0)> 



(22) 



is the overlap matrix. It should be noted that the overlap matrix is time- independent. 

The expectation values of the SGA observables are calculated in each realization (120]) of 
the weak measurement: 

M 

(j 4 ) = m)\ Ji !*(*)> = ■ (23) 

Averaging over the realizations gives the expectation values in the surrogate open-system 
dynamics (TjJ 



(24) 



Convergence of the simulation is checked on the level of Eq. (I24j) . 



2. Choice of numerical parameters 

A converged solution, Eq. (1241) . is independent of the size M of the computational basis, 
of the strength of the measurement 7 and of the size of the time-step dt in the Eq. (r2~TD . The 
convergence is obtained by increasing M and decreasing dt and 7. 
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a. Choice of the basis. A spin coherent state can be represented by a point on the 
phase-space sphere. It is expected that a state having high generalized purity can be repre- 
sented by a superposition of spin-coherent states localized in a small area of the phase-space. 
Solution of the sNLSE for an initial spin-coherent state has high generalized purity if the 
measurement is sufficiently strong 0,0]. Therefore, a localized computational basis is ex- 
pected to span the solution. Representation in terms of spin-coherent states is non-unique in 
general and the choice of the localized basis is non-unique, too. As a guideline for the choice 
of the basis we have used the following consideration. A set of iV + 1 spin-coherent states 
comprises the full basis of the Hilbert space, which can be represented by a distribution of 
N + 1 points on the sphere. A natural choice is a uniform distribution. Then a localized 
computational basis of M states can be mapped into M points uniformly distributed over 
an area A of the phase-space unit sphere, with A = 4nM/ (N+ 1). This mapping determines 
the parametrization of the computational basis. 

The value of M in the convergent solution of the sNLSE depends on the strength of 
measurement 7. It can be estimated balancing the rate of the derealization, induced by the 
Hamiltonian, and the rate of the localization, induced by the measurement. The derealiza- 
tion rate in the state prepared in the vicinity of the mean- field separatrix with U/ou — 2 is 
of the order of u 



271 ]. The localization rate in the state having generalized purity Psu(2)[ip], 
is of the order of 7X^=1 — (^) ^ ~ 7iV 2 (l — P su (2)[^]) |2j, which can be brought 

to the form ^MN using relation M < N (l — \J P B u(2)W\) derived in Ref.Q. Equating the 
rates we obtain 

M < — . (25) 

As 7 goes to zero M approaches its value in the corresponding unitary dynamics, which 
in the strong interaction regime may become of the order of N. The convergent solution of 
Eq.flU) is obtained at 7 < uj/(N\nN) [Sec. flHIB2cp ], The value of M in ([25]) is determined 
by 7 e = eu;/(iVm N), e 1. Therefore, 



b. Choice of the time-step. Eq. (fl5j) implies that the size of the time-step should 
be chosen to satisfy | |H ]■?/>) || dt~ l . Since adding a (time-dependent) phase to the 
state does not change the expectation values, the Hamiltonian can be substituted for 

15 



H — ( H } in the equations of motion. The condition of the size of the time-step becomes 

" \ /~\ 2 V /2 

H 2 } —(H) <C dt 1 , i.e., the size of the time-step must 



H - (H 



be smaller than the uncertainty of the Hamiltonian. The uncertainty of the Hamiltonian 
can be estimated in a localized state ip. The phase-space representation of the Hamiltonian 
quadratic in the operators Jj, as in Eq.flBJ), varies on the scale of unity, whereas the width 



of the phase-space representation of the localized state for M <C N is ~ a/ M/N < 1. As a 
consequence, the uncertainty of the Hamiltonian is of the order of ^H^ a/ M/N < ujy/MN 
and the time-step must satisfy 

uodt < (NM)~ 1/2 . (27) 

A similar estimation of other terms in Eg. (1151) . using condition 7 <C u/(Nha.N) for a 

convergent simulation [Sec. (1III B 2 cl) ] and Eqs. ([Hl) , leads to 

IniV , . 

Wdt ^ ~M~> (28) 

Since inequality ( j27j) is stronger than inequality ( l28l) . the scaling of the size of time-step of 
the simulation is determined by inequality (1271) . As a consequence, the time of computation 
scales as y/n, where n = N + 1 is the Hilbert space dimension of the system (Cf. Fig. [6]). 

c. Choice 0/7. In order that the effect of measurement on the dynamics of the target 
observables J r , J,, and J z can be neglected, it is sufficient that the contribution of the 
m e_t t0 I gr0 w th of t* 1. —ty Ell (()f ) - (5,) 2 ) ta the open- 
system evolution (jlj) is negligible on the time-scale of the simulation, compared to the 
uncertainty of the initial spin-coherent state. 

To estimate the contribution of the measurement to the growth of total uncertainty we 
consider dynamics of the single-particle observables under the action of measurement alone 

Ut)) = (U0))e- 2 ^. (29) 



The corresponding evolution of the total uncertainty is 

A( t ) = E((j. 2 )-(5.) 2 ) = ^(T + i)-^ W(0)>2e " 

j =1 \ / \ / i=1 



■4rft 
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The uncertainty in a spin-coherent state is minimal: A m « n = N/2. The effect of measurement 
can be neglected on a time-interval r if 



A(t) - A min (31) 

From Eqs.flUD and §T§ 



2 V2 / 2 2 4 v ; 2 

7riV 2 < — =3> 

7 « (2ATT)- 1 , (32) 

where assumption 1 <C A" is used. For the simulation on the time interval of order a; -1 In iV 
condition (1321) becomes 

Therefore, the open-system dynamics of the target observables converges to the original 
unitary evolution at the strength of measurement 7 e = c N ^ N , e < 1. 



3. Numerical results 



As an application of the algorithm we have performed the simulation of the system of 
N = 2 ■ 10 4 and N = 2 ■ 10 6 atoms. Panels (a) and (b) in Fig. ([3]) display dynamics of 
the target observables for two different realizations of the weak measurement governed by 
sNLSE (1131) . Each solution is converged with respect to the size M of the time-dependent 
basis. The size of the basis is M = 60 and the value of 7 is 0.3u/(Nln N). The generalized 
purity ([TBI) stays very close to unity as a consequence of the localization. The evolution of 
the expectation values of target observables in the two realizations of the weak measurement 
are identical on the time interval < t < 0.15co> _1 In N, corresponding to high generalized 
purity of the exact solution [Cf. FigfJ]. On this time interval the mean-field solution [Cf. 
Figf5] is a perfect approximation. At t > 0.3u; _1 In N the two realizations start to diverge. 
The computation of each realization takes about 30 seconds on a PC. 

Panel (c) shows the dynamics of the target observables averaged over the two realizations 
and the corresponding generalized purity. An important observation is the decrease of the 
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FIG. 3: (Color online). Panels (a) and (b) display solutions of sNLSE f)13[) . corresponding to two 
different realizations of the weak measurement. The normalized expectation values of the target 
observables J' x = jjjj x (blue), J' y = jfJ y (red) and J' z = jjJ z (green), Eq.([2]) and the generalized 
purity (black), Eq. (|16p . are plotted as functions of time. The system of N = 2 ■ 10 4 atoms is 
prepared in a spin-coherent state corresponding to the condensate placed in the left well of the 
trap. The value of the on-site interaction is chosen U/u = 2. Time is measured in units of lo" 1 In N. 
The size of the time-dependent basis M = 60. Panel (c) displays the average of the two realizations 
and the corresponding generalized purity. 

generalized purity. It implies the loss of the generalized purity in the exact dynamics and 
gives an estimate of the derealization time-scale. We conclude that some features of a 
many-body dynamics can be inferred from a small number of realizations. 

Averaging over 10 4 stochastic trajectories leads to a relative error of ~ 1% in the expec- 
tation values of the observables. Figure (j3J) presents the averaged results. The averaging has 
been performed over solutions of the sNLSE, corresponding to e = 7 iV In iV/u; in the range 
0.1 < e < 35. The averages are seen to converge as e decreases. The evolution on the time 
interval t < 0.3u; _1 lnA^ corresponds to the mean-field behavior. At t > 0.5u; _1 In TV the 
generalized purity decreases. Results of the simulation are compared to the exponentially 

n 

convergent nuenerica! solution 129], obtained by the sparse-matrix Cbebysbev propagation 
method [3|. The correspondence of the converged results of the two algorithms stays within 
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time in units of co~ 1 ln(N) 

FIG. 4: (Color online). N = 2 • 10 4 . Expectation values of the the target observables J' x = 
jjj x (blue solid), J' y = j^3 y (red solid) and J' 2 = jjJ z (green solid), Eq.([2]), averaged over 10 4 
realizations of the sNLSE CC| for the set of values {0.13, 0.26, 0.53, 1.06, 2.11, 4.23, 8.45, 16.9, 33.8} 
of e = jNlnN/cu. The expectation values and the corresponding generalized purity (black), 
Eq. (|16p . are plotted as functions of time. The initial state of the system and the value of U/u> are 
the same as in Fig. (3). Time is measured in units of w _1 ln(A r ). The results of the exponentially 
convergent computation using sparse-matrix Chebyshev techniques are presented for comparison 
(black dashed). 

1 — 2%. The size of the basis for the convergent solution on the interval < t < uj^ 1 InN is 
M s nlse — 60 and on the interval < t < l.'buo^ 1 In iV is M s ^lse — 180. The number of the 
time-steps for the propagation on the interval < t < 1.5c<j _1 In iV is 1.5 • 10 5 . 

FigJSl presents results of simulating the system of N = 2 • 10 6 atoms. The simulation 
is converged by averaging over 5 • 10 3 solutions of the sNLSE to ~ 1.5%. The size of the 
computational basis is M = 60. The number of timesteps is 8 ■ 10 5 . 

The scaling of the number of time-steps with the the Hilbert space-dimension n = N + 1 
was checked by finding the minimal number of time-steps k per unit time for fixed M, where 
the computation is still stable. The log-log plot for M = 60 in Fig. [6] shows the yjn scaling 
of k. Convergence of the simulations is obtained for the number of time-steps about 5 — 10 
times k per unit time. For N = 2 ■ 10 4 the simulation using the Chebyshev method is faster 
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FIG. 5: (Color online). N = 2 ■ 10 6 . Expectation values of the the target observables J' x = 
(blue), J'y = jfJy (red) and J' z = jjJ z (green) and the corresponding generalized purity (black) as 
functions of time. The initial state of the system and the value of U/uj are the same as in Fig. (3). 
The simulation is converged to ~ 1.5%, corresponding to averaging over 5000 realizations of the 
sNLSE (|13p . The size of the computational basis M = 90. 



due to the necessary averaging over the measurement realizations in the our method. For 
N = 2 • 10 6 our method is by orders of magnitude more efficient. 



IV. CONCLUSIONS 



The details and the first application of an algorithm for simulation of a many-body 
dynamics, generated by a su(2)-Hamiltonian, are reported. The algorithm implements idea, 
proposed in Ref.[2j, of using a surrogate open-system dynamics for simulation of a unitary 
evolution of the quantum many-body systems. The algorithm is applied to simulation of 
dynamics of iV = 2 • 10 4 and iV = 2 • 10 6 cold atoms in the double- well trap. The dynamics 
reflects a competition between the hopping rate of the atoms from well to well and the 
two-body repulsive interaction between the particles. The single-particle observables of the 
system are simulated and the convergence of the simulation is checked. 

The simulation is based on numerical solution of the stochastic Nonlinear Schrodinger 
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FIG. 6: (Color online). Empirical scaling of the minimal number of time steps k per unit time 
necessary for convergent simulation, as a function of the Hilbert space dimension n = N + 1. 
The estimate is performed by decreasing the number of steps until the computation becomes 
unstable. For convergent simulation it is sufficient to take 5 to 10 times k number of steps. The 
results correspond to the size of the computational basis M = 60. Stars and circles correspond 
to two different methods of performing time discretization. The dashed line (red) is a fit to 
ln(fc) = 0.51n(n) + c. 

Equation (sNLSE), which governs a single realization of the surrogate open-system dynamics. 
The open-system dynamics can be interpreted as a weak measurement of the elements of 
spectrum-generating su(2) algebra of the single-particle observables. It is shown that a 
range of the measurement strength exists, where on one hand the dynamics of the target 
observables is negligibly affected by the measurement but on the other hand the measurement 
is strong enough to localize solutions of the sNLSE in the associated phase-space. 

The localization leads to a drastic reduction of the computational complexity of the 
sNLSE compared to the original Schrodinger equation. The algorithm employs a time- 
dependent basis of the spin-coherent states. Asymptotically the number of states in the 
basis scales logarithmically with the Hilbert space dimension n = N + 1 for the simulation 
on the time-interval of order Inn. Therefore, the memory resources and the cost of an ele- 
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mentary time-step of the simulation are reduced by the factor ^ compared to conventional 
sparse-matrix algorithms for solving the corresponding Schrodinger Equation. The number 
of the time-steps per unit time for the convergent simulation on the fixed time interval scales 
as n 1 / 2 , i.e., is smaller by the factor n 1 / 2 than in conventional methods. This leads to the 
overall scaling of 0(\/n\n 3 n) of the computational resources for solving the sNLSE. The 
estimate of Ofn 1 / 2 ) scaling of the number of time-steps is rooted in the present linear imple- 
mentation of the algorithm. It is possible that the more flexible nonlinear implementation 
will lead to a better scaling. Solution of the sNLSE corresponds to a single realization of the 
surrogate open system dynamics. To obtain a reasonable accuracy of the simulation solu- 
tions of the sNLSE are averaged over ~ 10 3 realizations. Though in practice the averaging 
increases substantially the computational cost of the simulation, the number of realizations 
is independent on n and does not contribute to the asymptotic scaling. An important ad- 
ditional feature of the algorithm is that realizations of the open system dynamics can be 
simulated in parallel, which is advantageous from the computational standpoint. 

Simulating dynamics on the time-interval of order Inn is necessary to observe the gen- 
eralized purity loss, corresponding to depletion of the condensate, which is not seen in a 
mean-field approximation. Therefore, the Inn factor cannot be removed from the overall 
scaling of the algorithm and an efficient simulation of the depletion of the condensate is im- 
possible in the sense of our definition. It is conjectured that generally a logarithmic factor 
is present in the scaling of the computational resources with the Hilbert space dimension 
in a simulation of quantum many-body effects, i.e., effects that cannot be accounted for by 
a mean-field dynamics. This is because generically the mean-field approximation is broken 
when the generalized coherent state with respect to the SGA approaches the hyperbolic 
fixed point of the associated mean-field dynamics. The breakdown sets in when the distance 
to the fixed point becomes of the order of the total uncertainty of the generalized coherent 
state. The corresponding time-scale is of the order of the logarithm of the total uncertainty, 
which is of the order of the maximal weight of the Hilbert space representation of the SGA, 
increasing with the Hilbert space dimension. 

To summarize, the algorithm described in the present work reduces the memory resources, 
the cost of an elementary time-step of the simulation and the CPU time resources compared 
to conventional sparse-matrix approaches. As a consequence, asymptotically, it outperforms 

3/2 

conventional sparse-matrix computations by the factor ^— . It is conjectured that imple- 
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mentation of the algorithm to simulate many-body dynamics with other SGA will prove 
advantageous as well. 
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